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Abstract 

The network flow optimization approach is offered for Bayesian 
segmentation of gray-scale and color images. It is supposed image 
pixels are characterized by a feature function taking finite number of 
arbitrary rational values (it can be either intensity values or other 
characteristics of images) . The clusters of homogeneous pixels are de- 
scribed by labels with values in another set of rational numbers. They 
are assumed to be dependent and distributed according to either the 
exponential or the Gaussian Gibbs law. Instead traditionally used lo- 
cal neighborhoods of nearest pixels the completely connected graph of 
dependence of all pixels is employed for the Gibbs prior distributions. 

The methods developed reduce the problem of segmentation to 
the problem of determination of the minimum cut of an appropriate 
network. 
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1 Introduction 



The segmentation of images and restoration of degraded images are branches 
of image processing that are now extensively studied for their evident practi- 
cal importance as well as theoretical interest. There are many approaches to 
solution of these problems. We consider here methods of the Bayesian seg- 
mentation of images, which generalize methods of Bayesian image estimation 

|2|, H . Because of use of a prior information the methods of Bayesian image 
estimation would be methods of the first choice for many practical problems 
but unfortunately they are often difficult to compute. Until recently approxi- 
mations of the Gibbs estimators have been usually available rather then their 
exact values. For last decade the significant progress in the Gibbs estimation 
have been achieved. In particular, the methods of discrete optimization for 
the high- dimensional Gibbs estimation and segmentation have been obtained 
|| ||, ^, H, U- The methods developed allowed not only an efficient evaluation 
of the Gibbs estimates and but also their fast exact determination. 

In this paper the problem of image segmentation is considered as a prob- 
lem of cluster-analysis for the case of Gibbs prior distributions of clusters. 
There are enormous number of papers devoted to the problems of Bayesian 
cluster-analysis (see (II], [TTL [T2], |HJ). In some of them methods of classifica- 
tion with dependent clusters are presented Jl2|, [13| . We consider the problem 
of classification of observations yi,y2, ■■■ ,y n (f° r instance, image intensi- 
ties or texture characteristics of images etc.) with the feature function /(yj) 
that takes finite number of rational values. The clusters are identified by 
appropriate rational labels (note at once that in the presented below models 
the case of rationaly valued feature functions and cluster labels is reduced 
to the case of functions and labels taking only integer values). The clusters 
are supposed dependent and distributed according to the Gibbs field (such 
type models occur frequently in image processing). The labeling Gibbs field 
is specified on the directed fully connected graph of all pixels (usually only 
graphs of connected nearest neighbors pixels were considered. Often it was 
finite d-dimensional regular lattices). 

The discrete optimization methods that enable efficient determination of 
an exact solution of the segmentation problem are described. In the offered 
methods the problem of identification of the Gibbs classifiers is reduced, first, 
to the integer optimization problem and then to the problem of findind of 
the minimum cuts of special networks. The minimum cuts of the networks 
can be found by fast methods || [TJj that take into consideration the specific 



character of the networks and allow execution in concurrent mode. 
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2 Description of Models 



It can be easily seen that in models presented below the case of rationally 
valued feature functions and cluster labels is reduced to the case of inte- 
ger classification. So, let an image y = (yi,y 2 , ■ ■ ■ ,y n ) be with the feature 
function / that takes finite number integer values in Zl = {0, 1, . . . , L} and 
let fi = f( yi ). Let M = {to(0),to(1),... ,m{k)}, (0 < m(0) < m{\) < 
. . . < m(k) < L) be a set of allowable cluster labels. Suppose the image 
y is partitioned into k + 1 clusters (the number k + 1 instead of k is taken 
only to simplify formulas) and each cluster is specified by integer number 
j, (0 < j < k), as well as by an appropriate fixed integer label nij G M. 

The feature vector f = (/i, f 2 , . . . , f n ) of image y is considered as a 
random variable with the non-identical exponential 

P™exp(fi) = °i ■ eX P {-Mfi -ml} 

or with the non-identical Gaussian 

P™gausUi) = Ci ■ exp {-\(fi - m,) 2 } 

Gibbs distribution (A« > 0). 

The labels of clusters are supposed to be dependent random vari- 
ables distributed according to the Gibbs field. The segmented image x = 
{xi, x 2 , ■ ■ ■ , x n } is specified on the fully connected directed graph G = (V, E) 
with the set of vertices V = {0,1,... ,n} and the set of directed arcs 
E = | i,j £ V}. It takes values in the space of labels M v , i.e. 

Xi G A4, i G V, and is either the exponential or the gaussian form 

p exp (m) = c-exp{- ^ Pijlm- rrijl}, 
p gaus {m) = c ■ exp{- ^ Pij{ m i ~ m j) 2 }i 

where c is the norming quantities (here and below different constant will be 
denoted by the same letter), the vector m = {mi, m 2 , . . . , m n } and param- 
eters Pij > 0. 

Remark 1. (i) The Gibbs models considered are extentions of the Ising 
model. 

(ii) Note, that for some image processing problems the prior distribution 
p exp {x) is more preferable than Gaussian Gibbs one because of smaller blur- 
ring effect. Moreover, its identification turns out far more computationally 
efficient. 
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For the model with the exponential conditional distribution p™ exp {fi) and 
the exponential prior p exp (m) the Gibbs classifier is equal to 



m exp = arg min < A; \fr - rrii\ + Ajl^i ~ m j\ \ , (!) 

and respectively, for the model with the Gaussian conditional distribution 
Pi^gaus(fi) and the Gaussian prior p gaU s(m) the Gibbs classifier is of the form 



m g aus = arg min < V] Ai(/ f - rrii) 2 + V" Pi,j(mi - mj) 2 > . (2) 

m£M v * — ' z — ' ' 



In spite of clear posing the problem of finding exact values of the Gibbs 
classifiers m exp and m gaus for large samples is rather complicated. So, for in- 
stance, computer images, which are frequent objects of classification, usually 
consist up to 2 18 variables or more. Nevertheless, it turned out possible to 
compute efficiently (during polynomial time on the sample size n and number 
of clusters k) both of these classifiers. Theoretically, run-time for the first 
classifier does not exceed ckn 3 and for the second one is less than c(kn) 3 . For 
many applied problems real run-time was even of order 0(nk). The efficient 
computation of the classifiers for the mixed model with the exponential con- 
ditional distributions p™ exp (fi) and the Gaussian prior p gaus (m) as well as for 
the mixed model with the Gaussian conditional distributions P™ gaus {fi) and 
the exponential prior p exp (m) is also available. 



3 Computation of the Optimal Classifiers 

In the case of two classes (k = 1) the Gibbs classifier m exp and m gaus coin- 
cide (since any Boolean variable b satisfies the identity b 2 = \b\). They can 
be evaluated by the network flow optimization methods 0, [15| . Moreover, 
in 1989 Greig, Porteous and Seheult developed a heuristic network flow 
algorithm that is especially efficient for estimating the Boolean Gibbs estima- 
tor. These authors posed also the problem for the case more than 2 clusters. 
Recently we have described the multiresolution network flow minimum cut 
algorithm |J that allows exact computation of Boolean classifiers as well as 
developed algorithms of computation of the mentioned Gibbs classifiers in 
general case. It turned out the network flow optimization methods can be 
used to identify m exp and m gaus even when k > 1. 
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Identification of rh 



exp 



Denote the function to be minimized by 

Ui{m) = ) j \j\fj - mj\ + Pijmi-mjl, (m G M v ). 

i&V (i,j)eE 

The idea of the method is to represent the vector of labels of clusters m by 
the integer valued linear combination of Boolean vectors and then reduce 
the problem of integer minimization of the function U\ (m) to the problem of 
Boolean minimization. The problem of Boolean minimization can be solved 
by the network flow optimization methods. 

Let for arbitrary integers \i and v the indicator function l(p> v ) be equal 
to 1 if n > v and be equal to otherwise. For any /i G M. and Boolean 
variables x{l) — 1^ > ^ such that x(l) > x{2) > . . . > x{k) the identity 



fx = m (0) + 5^(m(Z) - m(l - l))x(l) 



(3) 



i=i 



is valid, and vice versa, any non-increasing sequence of the Boolean variables 
x(l) > x{2) > . . . > x(k) specifies the label \i G M. by the formula (0). By 
analogy, the feature functions fi are represented as sums /, = J2t=i /«( r ) °^ 
non-increasing sequence of the Boolean variables > /i(2) > . . . > fi(L). 

Let for I = 1 -T- k the vector x(Z) = (xi(l) , X2(l) , ■ ■ ■ ,x n (l)) be Boolean, 
the vector z(7) = (zi(l), . . . , z n {l) be with coordinates 



Ziil) 



m{l) — m{l — 1) 



m(l) 



(i 



n 



r=m(i-l)+l 



and the norm of the vector |x| = Y^h=i \ x i\-> then the following proposition is 
satisfied. 

Proposition 1. For any integers v G M. and f) G the equality 

W - fi\ 



m(0) 

m(0) - J2 fi(T) 



T=l 



+ 



m(l) 

(m(l)-m(l-l))lr y ^ U 

V ' T=m{l-1)+1 



E 



Z=l 



T=m(fc)+1 



ZioZds true. Therefore, for any feature vector f G Z\ and the Boolean vector 



x(/) 



vl(mi > 0' ■'■(ma > 0' 



' J -(m„ > Z 



)) 
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the function Ux(m) can be written in the form 

k 

LMm) = E( m ^ " m ( Z " WWQ), ( 4 ) 

where for n- dimensional Boolean vector b functions 

u(l,b) = y^Ail^(/) - 6f| + Pijlk-bjl. 

Denote by 

x(Z) = axgminw(7, b), (Z = 1 4- (5) 

b 

Boolean solutions that minimize the functions u(l,h). For two vectors v 
and w we will write v > w if all corresponding pairs of their coordinates 
satisfy the inequality Vi > Wi, (i — 1 -f- n) and will write u ^ t/7 if there 
exist at least two different pairs such that v,i > uii and Vj < Wj. Note that 
z(l) > z(2) > . . . > z(k), and what is more, for some integer 1 < x < k 
their coordinates satisfy the following condition 

Zi {\) = Zi {2) = ... = *(*- 1) = 1, 1, 

Zi(x+1) = ... = Zi(k) = 0. 

It is easy to show that in general the solutions x(Z) of the (H) are not or- 
dered. Nevertheless, without fail there is at least one non-increasing sequence 
x(l) of solutions of (||). 

Theorem 2. There is a non-increasing sequence x(l) > x(2) > . . . > x(&) 
of solutions of (fjj. 

Some structural properties of the set of solutions x(7) are presented in 

Corollary 3. For integer 1 </'</"< k and the sequence of vectors z(l) > 
z(2) > . . . > z(k) the following properties are valid: 

(i) If x(/') is any solution of (j^), then there is a solution x(/") so that 
x(7') > x(/"), and vice versa, if x(7") is any solution of (fjj, then there is a 
solution x(T) so that x(Z') > x(7"). 

(ii) For each 1 < I < k the set of solutions {x(/)} has the minimal x(7) 
and the maximal x(/) elements. 

(Hi) The set of minimal and maximal elements are ordered, i.e. x(i) > 
. . . > x(/c) andx(l) > ... > x(A;). 
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Sentence (i) follows immediately from Theorem 0. Sentence (ii) is de- 
duced from (i) considered for I' = I", property (Hi) follows from (ii) and 
definition of the minimal and the maximal elements. 

For two Boolean vectors b' and b" let the vector b = b' A b", re- 
spectively, b = b' V b" be with coordinates = min{Z^,6"}, respectively, 
with coordinates \ = max{6-,Z>"}. If x(l), . . . , x(fc) is any unordered se- 
quence of solutions of (J5|) the ordered sequence of solutions can be derived 
from it by the logical operation A, V like one-dimensional variational series 
x (i) > x (2) > • • • > X(fe). But it is easy to see the sum of ordered solutions 
x(Z) is a solution of ([]]). 

Proposition 4. 7/x(l) > x(2) > . . . > x(/c) is a sequence of ordered solu- 
tions of (^) then the sum 

k 

m = m(0) + V}(ra(Z) - mil - l))x(Z) 

1=1 

minimizes Ui(m). 

The problem of computing Boolean solutions x(Z) is familiar in the dis- 
crete optimization j| . It is equivalent to identification of the minimum cuts 
for specially built networks. There are fast algorithms to compute them 
IH, W< 



Identification of m ffaus 

Now denote the function to be minimized by 

Z7 2 (m) = J2 X ^h ~ m ^ + Yl ~ m J') 2 ' ( m E MV }- 

To find a solution m gaus that minimizes the function ^(m) the represen- 
tation of the vector of cluster labels m by the integer valued linear combi- 
nation of Boolean vectors is used once more. Then the problem of integer 
minimization of the function ^(m) is reduced to the problem of Boolean 
minimization. 

Denote for brevity = /— m(0), (i G V), ai = m(l)—m(l—l), (Z = 1-rk). 
The vector m can be represented by the formula (^) as the linear combination 
m(0) + Ym=i a z x (0 °f Boolean vectors x(Z) = (xi(Z), x 2 (Z), ... ,x n (l)), and 
the function Z7 2 can be written in the form 

(k v 2 / fc \ 2 

gi - aixjjt) j + ^ ( q / (^(o ~ x j (o) ) • 
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Let di >T = aia T , (l,r = 1 k) and j3 i}i — 0, (i e V) , then ^(m) = 
X^ev + -P( x (l)> • • • ? x (^))i where the polynomial of Boolean variables 
P(x(l), . . . , x(fc)) after cancellation is written as 



P(x(l),...,x(fc)) = 

k r i 

^2 K&I ~ 2 ^i9i a i - ai{m(k) - m(0) - aj) ^(Aj + #,-,0 + 
r fc 

Note that it has the same points of minimum as [^(m). Let us consider 
another polynomial of Boolean variables 



Q(x(l),... ,x(fc)) = 

^2 S ~~ 2A ^ a ' ~~ ai(m(k) - m(0) - aj) ^(Aj + 
iav i=i L jev 

2 S A *+X^-+^) E ^(o + 



Xi(Z) + 



i€V 



KKKfc 



£ A,j ^af(x l (/)-x,(/)) 2 +^^ T [(^(/)-^(r)) 2 +(x,(/)-x J (r)) 2 ] 



such that Q(x(l), . . . , x(fc)) > P(x(l), . . . , x(fc)) and which differs from P 
by the term ^2 1<T<l<k di jT Xi(l) in the second line. 
Denote by 

(q*(l),q*(2),... ,q*(A;)) = argmin x(1))X(2)) ... !x(fe) Q(x(l), . . . ,x(fc)) 

any collection of Boolean vectors that minimizes Q(x(l), . . . , x(/c)). Without 
fail q*(l) > q*(2) > ... > q*(L — 1). This feature allows expressing solutions 
of the initial problem as m gaus = J2t=i 

Theorem 5. Any collection (q*(l), q*(2), . . . , q*(L — 1)) that minimizes Q 
forms the nonincreasing sequence. 

The polynomials P and Q have the same set of ordered solutions and, 
therefore, each solution m gaus is specified by the formula m gaus = Y^=i q*(0- 
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Theorem [5] allows determination of the classifier m gaus by the Boolean 
minimization of the polynomial Q. Unlike P this polynomial can be mini- 
mized directly by the minimum network cut algorithms 0, |S| . The appropri- 
ate network is described in ||. 

4 Applications 

Here we show several numerical tests as well as a resust of segmentation of 
a real 3D US-imade of size 196 x 215 x 301. 

The original 2D gray-scale image in Figure la (see next page) was cor- 
rupted by Gaussian random noise (Figure lb). The results of restoration by 
the extened Ising model are placed in Figure lc,ld and by the classical Ising 
model are depicted in Figure le,lf. 

In Figure 2a the slice of original 3D US-image of the thyroid gland is 
depicted. Its contour that was done by expert is drown in Figure 2a. The 
corresponding slice of 3D segmentation of the original image by the extended 
Ising model are placed in Figure 2c, d. The full segmentation of 3D 196 x 
215 x 301 image takes about 40min of processor Pentium-Ill 800. 
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(a) 



(b) 





(d) 



Figure 2 



5 Conclusion 

In the paper the Bayesian methods of segmentation and estimation of gray- 
scale and color images are presented. For both of them it is supposed feature 
function of images and labels of segments take finite number of rational 
(possibly, different) values and they are distributed according to either the 
exponential Gibbs or the Gaussian Gibbs distribution. The numerical tests 
showed the methods developed allow solution problems of practical segmen- 
tation and Gibbs estimation of images of large sizes. 
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